Este documento contém a sétima especificação do VAR que eu testei. Ela é a mesma variação do VAR 1, porém com diagnóstico de convergência.
Para as transformações que foram feitas nos dados e a motivação da escolha das variáveis, sugiro ver meu arquivo de descritivas.
Para fazer o BVar utilizei o pacote bvarsv do Fabian Kruger. Em particular, esses dois arquivos me ajudaram muito: bvarsv: An R implementation of the Primiceri (2005) model for macroeconomic time series e Replication of figures in Del Negro and Primiceri (2015).
Para fazer o diagnóstico da cadeia de Markov (algo que não estava presente nos outros VARs), usei o pacote coda junto com …
Se o mirror 10 (UFRJ) estiver fora do ar, tente mudar ind = 10 para ind = 1. Mesma coisa se o BETS não instalar.
chooseCRANmirror(graphics = FALSE, ind = 10)
if (!require("pacman")) install.packages("pacman")
pacman::p_load(ggplot2, forecast, BETS, seasonal, seasonalview, bvarsv, lubridate, zoo, stargazer, gridExtra, reshape2, ggfortify, RColorBrewer, scales, coda, foreach, doParallel, knitr, grid, ggpubr)
Algumas funções auxiliares:
## Peguei essa função pronta
matplot2 <- function(...){
matplot(..., type = 'l', lty = 1, lwd = 2, bty = "n", ylab = "")
}
# Essa função calcula os quantis e organiza 1 quantil, média, 3 quantil.
stat.helper <- function(z) c(mean(z), quantile(z, c(0.16, 0.84)))[c(2,1,3)]
# Paleta de cores
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
cols1 <- cbPalette[c(2,4,2)]
cols2 <- cbPalette[c(2,4,6)]
A primeira tentativa foi com dados de janeiro de 2003 a outubro de 2017 e é igual à especificação do primeiro VAR. Minha intenção era testar os diagnósticos e ver se o computador consegue salvar a séries completas do amostrados de Gibbs.
Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do trabalho (7620)Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do capital (7621)Passei um filtro para tirar sazonalidade, sem transformação prévia nos dados (ajuste default na função seas exceto no transformation que ficou igual a NONE).
Taxa de juros - Selic acumulada no mês (4390). Passei de taxa mensal para anual usando a fórmula: \(\left((1+tx/100)^12 -1\right)*100\) e passei um filtro para tirar sazonalidade, sem transformação prévia nos dados (ajuste default na função seas exceto no transformation que ficou igual a NONE).Índice nacional de preços ao consumidor-amplo (IPCA) (433). Como ele está em variação percentual mensal, foi necessário transformar para acumulado dos últimos 12 meses utilizando a fórmula: \[IPCA_i = \left[\left(\prod\limits_{j=i-11}^i \left(\frac{IPCA_j}{100}+1\right) \right) -1\right]*100\]Índice de Atividade Econômica do Banco Central (IBC-Br) - com ajuste sazonal (24364). Fiz a diferença da série em log, igual o câmbio.Taxa de câmbio - Livre - Dólar americano (venda) - Fim de período - mensal (3696). Como o Mumtaz disse que usou log diferença, eu calculei o log e tirei primeira diferença nessa série.#rm(list = ls())
# Auxiliary variables, so I don't need to bother when something changes
inicio <- "2003-02-01"
fim <- "2017-10-31"
# Don't mess with this code
inicio_cambio <- paste(seq(as.Date(inicio), length = 2, by = "-1 month")[2]) # 1 mês antes do início das outras séries
inicio_ipca <- paste(seq(as.Date(inicio), length = 12, by = "-1 month")[12]) # 12 meses antes do início das outras séries
trabalho <- BETS.get("7620", from = inicio, to = fim)
capital <- BETS.get("7621", from = inicio, to = fim)
capital_trabalho <- capital/trabalho
# Usando ndiffs(ibcbr,test="adf",alpha = 0.1) se encontra que o ibcbr é não estacionário.
# Como todas as outras séries são estacionárias, eu vou fazer também a diferença do log, igual no câmbio
ibcbr_raw <- BETS.get("24364", from = inicio_cambio, to = fim)
ibcbr <- diff(log(ibcbr_raw), 1)
cambio_raw <- BETS.get("3696", from = inicio_cambio, to = fim)
cambio <- diff(log(cambio_raw), 1)
selic_4390 <- BETS.get("4390", from = inicio, to = fim)
selic <- ((1+selic_4390/100)^(12)-1)*100
ipca_raw <- BETS.get("433", from = inicio_ipca, to = fim)
ipca_acum <- ipca_raw/100 + 1
ipca <- vector()
final <- length(ipca_acum)
for (i in 12:final){
ipca[(i-11)] <- (prod(ipca_acum[(i-11):i])-1)*100
}
ano <- as.numeric(substr(inicio, start = 1, stop = 4))
mes <- as.numeric(substr(inicio, start = 6, stop = 7))
dia <- as.numeric(substr(inicio, start = 9, stop = 10))
ipca <- ts(ipca, start = c(ano, mes, dia), frequency = 12) # Date format YYYY MM DD
m <- seas(x = capital_trabalho)
capital_trabalho2 <- final(m)
m <- seas(x = selic, transform.function = "none")
selic2 <- final(m)
df1 <- data.frame(seq(as.Date(inicio), length = length(ipca), by = "1 month"), capital_trabalho2, selic2, ibcbr, cambio, ipca)
names(df1) <- c("Data", "Capital_trabalho", "Selic", "IBCBr", "Cambio", "IPCA")
df2 <- melt(data = df1, id.vars = "Data")
cores <- brewer.pal(5, "Dark2")
# Gráficos individuais
p1 <- ggplot(df2[which(df2$variable == "Capital_trabalho"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[1])+
scale_y_continuous(name="Capital trabalho") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p1 <- p1 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p2 <- ggplot(df2[which(df2$variable == "Selic"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[2])+
scale_y_continuous(name="Selic (%a.a.)") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p2 <- p2 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p3 <- ggplot(df2[which(df2$variable == "IBCBr"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[3])+
scale_y_continuous(name="IBC-Br") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p3 <- p3 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p4 <- ggplot(df2[which(df2$variable == "Cambio"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[4])+
scale_y_continuous(name="Tx. Câmbio (var)") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p4 <- p4 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p5 <- ggplot(df2[which(df2$variable == "IPCA"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[5])+
scale_y_continuous(name="IPCA (acum. 12m.)") +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
theme_bw()
p5 <- p5 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
grid.arrange(p1, p2, p3, p4, p5, ncol=1, nrow = 5)
rm(capital, capital_trabalho, trabalho, ipca_acum, ipca_raw, ibcbr_raw, selic_4390, selic, cambio_raw)
rm(p1, p2, p3, p4, p5, df2)
### Descriptives
descriptives <- matrix(NA, nrow = 8, ncol = (ncol(df1)-1))
rownames(descriptives) <- c("Observações", "Mínimo", "1o quartil",
"Média", "Mediana", "3o quartil", "Máximo",
"Desv. Pad.")
colnames(descriptives) <- names(df1)[-1]
desc <- function(x) {
n <- length(x)
minimum <- min(x, na.rm = TRUE)
first_q <- quantile(x, 0.25, na.rm = TRUE)
media <- mean(x, na.rm = TRUE)
mediana <- median(x, na.rm = TRUE)
third_q <- quantile(x, 0.75, na.rm = TRUE)
maximum <- max(x, na.rm = TRUE)
std <- sd(x, na.rm = TRUE)
return(list(n = n, minimum = minimum, first_quar = first_q, media = media, mediana = mediana, third_quar = third_q, maximum = maximum, std = std))
}
for (i in 1:8){
descriptives[i, 1] <- round(as.numeric(desc(df1[,2])[i]),4)
descriptives[i, 2] <- round(as.numeric(desc(df1[,3])[i]),4)
descriptives[i, 3] <- round(as.numeric(desc(df1[,4])[i]),4)
descriptives[i, 4] <- round(as.numeric(desc(df1[,5])[i]),4)
descriptives[i, 5] <- round(as.numeric(desc(df1[,6])[i]),4)
}
descriptives[1,] <- as.integer(descriptives[1,])
descriptives <- data.frame(descriptives)
stargazer(descriptives, summary=FALSE, header = TRUE, type = 'html')
| Capital_trabalho | Selic | IBCBr | Cambio | IPCA | |
| Observações | 177 | 177 | 177 | 177 | 177 |
| Mínimo | 0.190 | 6.707 | -0.031 | -0.149 | 2.456 |
| 1o quartil | 0.436 | 10.305 | -0.003 | -0.027 | 4.832 |
| Média | 0.526 | 12.953 | 0.002 | -0.0004 | 6.497 |
| Mediana | 0.493 | 12.263 | 0.002 | -0.005 | 6.037 |
| 3o quartil | 0.570 | 14.416 | 0.007 | 0.022 | 7.138 |
| Máximo | 1.305 | 26.608 | 0.019 | 0.158 | 17.235 |
| Desv. Pad. | 0.163 | 4.105 | 0.008 | 0.045 | 2.799 |
#stargazer(descriptives, summary=FALSE, header = TRUE, type = 'latex')
rm(descriptives, df1)
Inicialmente ajustei um modelo com uma defasagem e que utiliza as primeiras 24 observações para fazer a estimativa de MQO para jogar na priori.
var1 <- cbind(capital_trabalho2, selic2, ibcbr, cambio, ipca) # The bvar function does not allows data.frames
set.seed(1)
nburn. <- 5000
nrep. <- 50000
fit1 <- bvar.sv.tvp(var1, p=1, tau = 24, nburn = nburn., nrep = nrep., nf = 1, thinfac = 1)
## [1] "2018-01-26 06:55:45 -- now starting MCMC"
## [1] "2018-01-26 06:56:51 -- now at iteration 5000"
## [1] "2018-01-26 06:59:11 -- now at iteration 15000"
## [1] "2018-01-26 07:01:26 -- now at iteration 25000"
## [1] "2018-01-26 07:03:43 -- now at iteration 35000"
## [1] "2018-01-26 07:06:00 -- now at iteration 45000"
## [1] "2018-01-26 07:08:17 -- now at iteration 55000"
## Colocando as datas como strings
tm1 <- as.yearmon(time(var1))
# Eixo x
xax <- time(var1)
# Perde uma observação pelo lag e outras 24 pela amostra para estimar por MQO, então começa no tempo 26
xax <- xax[26:177]
# Marcas verticais
gp <- seq(1999, 2017, 1)
Peguei a explicação daqui daqui: A FIR estima o impacto de um choque unitário em algum elemento de \(\varepsilon_t\). Os elementos de \(c_t\), \(B_{j,t}\), \(A_t\) e \(\Sigma_t\) em (1) são mantidos fixos em seus valores em \(t\). Note que a ordem das variáveis importa e deve ser justificada por argumentos econômicos. As FIR contém os percentis 5, 25, 50, 75 e 95.
# Função Impulso resposta da Selic na razão capital-trabalho
ira <- impulse.responses(fit1, impulse.variable = 2, response.variable = 1)
# Função Impulso resposta da Selic no IBC-Br
ira <- impulse.responses(fit1, impulse.variable = 2, response.variable = 3)
# Função Impulso resposta da Selic no câmbio
ira <- impulse.responses(fit1, impulse.variable = 2, response.variable = 4)
# Função Impulso resposta Selic no IPCA
ira <- impulse.responses(fit1, impulse.variable = 2, response.variable = 5)
# Função Impulso resposta do IBC-Br na razão capital trabalho
ira <- impulse.responses(fit1, impulse.variable = 3, response.variable = 1)
# Função Impulso resposta do Cambio na razão capital trabalho
ira <- impulse.responses(fit1, impulse.variable = 4, response.variable = 1)
# Função Impulso resposta do IPCA na razão capital trabalho
ira <- impulse.responses(fit1, impulse.variable = 5, response.variable = 1)
rm(var1, ira)
# Saving parameters
beta_ct <- parameter.draws(fit1, type = "lag1", row = 1, col = 1)
beta_selic <- parameter.draws(fit1, type = "lag1", row = 2, col = 2)
beta_ibcbr <- parameter.draws(fit1, type = "lag1", row = 3, col = 3)
beta_cambio <- parameter.draws(fit1, type = "lag1", row = 4, col = 4)
beta_ipca <- parameter.draws(fit1, type = "lag1", row = 5, col = 5)
sd_ct <- parameter.draws(fit1, type = "vcv", row = 1, col = 1)
sd_selic <- parameter.draws(fit1, type = "vcv", row = 2, col = 2)
sd_ibcbr <- parameter.draws(fit1, type = "vcv", row = 3, col = 3)
sd_cambio <- parameter.draws(fit1, type = "vcv", row = 4, col = 4)
sd_ipca <- parameter.draws(fit1, type = "vcv", row = 5, col = 5)
# salva o fit em um csv e limpa da memória
saveRDS(fit1, "C:\\Users\\aisha\\Documents\\Var da Aisha\\Var7\\fit1.rds")
#fit1 <- readRDS("C:\\Users\\aisha\\Documents\\Var da Aisha\\Var7\\fit1.rds")
rm(fit1)
# Convert stuff into mcmc objects
mcmc_beta_ct <- mcmc(beta_ct, start = 1, end = 50000)
mcmc_beta_selic <- mcmc(beta_selic, start = 1, end = 50000)
mcmc_beta_ibcbr <- mcmc(beta_ibcbr, start = 1, end = 50000)
mcmc_beta_cambio <- mcmc(beta_cambio, start = 1, end = 50000)
mcmc_beta_ipca <- mcmc(beta_ipca, start = 1, end = 50000)
mcmc_sd_ct <- mcmc(sd_ct, start = 1, end = 50000)
mcmc_sd_selic <- mcmc(sd_selic, start = 1, end = 50000)
mcmc_sd_ibcbr <- mcmc(sd_ibcbr, start = 1, end = 50000)
mcmc_sd_cambio <- mcmc(sd_cambio, start = 1, end = 50000)
mcmc_sd_ipca <- mcmc(sd_ipca, start = 1, end = 50000)
Os gráficos abaixo contém os betas (esq) e volatilidade (dir) ao longo do tempo, com o 16o e o 68o percentis (como foi feito no artigo do Primiceri).
# Gráfico do desvio padrão do resíduo da razão capital trabalho
x1 <- t(apply(sqrt(sd_ct), 2, stat.helper))
x1betact <- t(apply(beta_ct, 2, stat.helper))
#x1beta <- cbind(desc_betact$quantiles[,1], desc_betact$statistics[,1], desc_betact$quantiles[,5]) # ficou muito amplo mas meu coração não deixou eu apagar essa linha
desc_betact <- summary(mcmc_beta_ct)
desc_sdct <- summary(mcmc_sd_ct)
old.par <- par(mfrow=c(1,2))
matplot2(x = xax, y = x1betact, ylim = c(min(desc_betact$quantiles[,1]), max(desc_betact$quantiles[,5])), col = cols1, main = "Betas da Razão Capital/Trabalho", xlab = "Tempo")
matplot2(x = xax, y = x1, ylim = c(sqrt(min(desc_sdct$quantiles[,1])), sqrt(max(desc_sdct$quantiles[,5]))), col = cols1, main = "Volatilidade da Razão Capital/Trabalho", xlab = "Tempo")
#matplot2(x = xax, y = x1, ylim = c(0.04, 0.065), col = cols1, main = "Volatilidade da Razão Capital/Trabalho", xlab = "Tempo")
#matplot2(x = xax, y = x1, ylim = c(min(x1[,1])-2*sd(x1[,1]), max(x1[,3])+2*sd(x1[,3])), col = cols1, main = "Volatilidade da Razão Capital/Trabalho", xlab = "Tempo")
par(old.par) # resets par()
rm(desc_betact, desc_sdct)
# Gráfico do desvio padrão do resíduo da selic
x2 <- t(apply(sqrt(sd_selic), 2, stat.helper))
x2betaselic <- t(apply(beta_selic, 2, stat.helper))
desc_betaselic <- summary(mcmc_beta_selic)
desc_sdselic <- summary(mcmc_sd_selic)
old.par <- par(mfrow=c(1,2))
matplot2(x = xax, y = x2betaselic, ylim = c(min(desc_betaselic$quantiles[,1]), max(desc_betaselic$quantiles[,5])), col = cols1, main = "Betas da Selic", xlab = "Tempo")
matplot2(x = xax, y = x2, ylim = c(sqrt(min(desc_sdselic$quantiles[,1])), sqrt(max(desc_sdselic$quantiles[,5]))), col = cols1, main = "Volatilidade da Selic", xlab = "Tempo")
par(old.par) # resets par()
rm(desc_betaselic, desc_sdselic)
# Gráfico do desvio padrão do resíduo do ibcbr
x3 <- t(apply(sqrt(sd_ibcbr), 2, stat.helper))
x3betaibcbr <- t(apply(beta_ibcbr, 2, stat.helper))
desc_betaibcbr <- summary(mcmc_beta_ibcbr)
desc_sdibcbr <- summary(mcmc_sd_ibcbr)
old.par <- par(mfrow=c(1,2))
matplot2(x = xax, y = x3betaibcbr, ylim = c(min(desc_betaibcbr$quantiles[,1]), max(desc_betaibcbr$quantiles[,5])), col = cols1, main = "Betas do IBC-Br", xlab = "Tempo")
matplot2(x = xax, y = x3, ylim = c(sqrt(min(desc_sdibcbr$quantiles[,1])), sqrt(max(desc_sdibcbr$quantiles[,5]))), col = cols1, main = "Volatilidade do IBC-Br", xlab = "Tempo")
par(old.par) # resets par()
rm(desc_betaibcbr, desc_sdibcbr)
# Gráfico do desvio padrão do resíduo do cambio
x4 <- t(apply(sqrt(sd_cambio), 2, stat.helper))
x4betacambio <- t(apply(beta_cambio, 2, stat.helper))
desc_betacambio <- summary(mcmc_beta_cambio)
desc_sdcambio <- summary(mcmc_sd_cambio)
old.par <- par(mfrow=c(1,2))
matplot2(x = xax, y = x4betacambio, ylim = c(min(desc_betacambio$quantiles[,1]), max(desc_betacambio$quantiles[,5])), col = cols1, main = "Betas da Variação da Tx. de Câmbio", xlab = "Tempo")
matplot2(x = xax, y = x4, ylim = c(sqrt(min(desc_sdcambio$quantiles[,1])), sqrt(max(desc_sdcambio$quantiles[,5]))), col = cols1, main = "Volatilidade da Variação da Tx. de Câmbio", xlab = "Tempo")
par(old.par) # resets par()
rm(desc_betacambio, desc_sdcambio)
# Gráfico do desvio padrão do resíduo do ipca
x5 <- t(apply(sqrt(sd_ipca), 2, stat.helper))
x5betaipca <- t(apply(beta_ipca, 2, stat.helper))
desc_betaipca <- summary(mcmc_beta_ipca)
desc_sdipca <- summary(mcmc_sd_ipca)
old.par <- par(mfrow=c(1,2))
matplot2(x = xax, y = x5betaipca, ylim = c(min(desc_betaipca$quantiles[,1]), max(desc_betaipca$quantiles[,5])), col = cols1, main = "Betas do IPCA acum. 12 meses", xlab = "Tempo")
matplot2(x = xax, y = x5, ylim = c(sqrt(min(desc_sdipca$quantiles[,1])), sqrt(max(desc_sdipca$quantiles[,5]))), col = cols1, main = "Vol. do IPCA acum. 12 meses", xlab = "Tempo")
par(old.par)
rm(desc_betaipca, desc_sdipca)
Estou me baseando nos seguintes linkis: * www.patricklam.org/teaching/convergence_print.pdf para traceplot, density plot, running mean plot (é o plot da iteração contra a média das observações até a iteração (ele está indo de 10 em 10 para frente na média)) * https://theoreticalecology.wordpress.com/2011/12/09/mcmc-chain-analysis-and-convergence-diagnostics-with-coda-in-r/
Um dos diagnósticos é plotar a cadeia (usando plot(mcmc_object)). Ele retorna o gráfico dos valores (um gráfico de linha - quanto mais ruidoso melhor) e um traceplot que é o gráfico da distribuição. Se eu não estou enganada, os meus Betas tem distribuição à posteriori normal e falta ver a distribuição dos sigmas. Se a densidade das matrizes é inversa wishart eles tem distribuição gamma, mas acho que não tem fórmula fechada para \(\Sigma\) - preciso ver no artigo do Primeri depois.
Outro diagnóstico das médias envolve fazer o gráfico da média usando uma janela móvel para ver se a posteriori cai em alguma região e fica empacada.
Como eu tenho muitos instantes de tempo (e portanto uma quantidade enorme de betas e de sigmas), eu selecionei arbitrariamente 15 betas (que é a divisão inteira de t por 10 - então estou plotando 10% dos parâmetros para cada \(y\)).
sequencia <- seq(10, dim(beta_ct)[2], 10)
sequencia2 <- seq(1:dim(beta_ct)[1])
#####################################
# Para a razão capital trabalho
#####################################
#############
# Betas
#############
# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
plot(mcmc_beta_ct[,i], auto.layout = T, main = paste("Beta em t=", as.character(i), "da Razão Cap./Trab."))
}
par(old.par)
# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
df1 <- data.frame(sequencia2,mcmc_beta_ct[,i])
names(df1) <- c("Iteracao", "Beta")
ggplot(df1, aes(Iteracao, Beta)) +
geom_line(aes(y = rollmean(Beta, 10, na.pad = TRUE))) +
scale_y_continuous(name = paste("Média móvel (10 valores) de Beta para t=", as.character(i))) +
theme_bw()
}
par(old.par)
##############
# Volatilidade
##############
# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
plot(mcmc_sd_ct[,i], auto.layout = T, main = paste("D.P. em t=", as.character(i), "da Razão Cap./Trab."))
}
par(old.par)
# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
df1 <- data.frame(sequencia2,sqrt(mcmc_sd_ct[,i]))
names(df1) <- c("Iteracao", "DP")
ggplot(df1, aes(Iteracao, DP)) +
geom_line(aes(y = rollmean(DP, 10, na.pad = TRUE))) +
scale_y_continuous(name = paste("Média móvel (10 valores) do desvio padrão para t=", as.character(i))) +
theme_bw()
}
par(old.par)
#####################################
# Para a Selic
#####################################
#############
# Betas
#############
# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
plot(mcmc_beta_selic[,i], auto.layout = T, main = paste("Beta em t=", as.character(i), "da Selic"))
}
par(old.par)
# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
df1 <- data.frame(sequencia2,mcmc_beta_selic[,i])
names(df1) <- c("Iteracao", "Beta")
ggplot(df1, aes(Iteracao, Beta)) +
geom_line(aes(y = rollmean(Beta, 10, na.pad = TRUE))) +
scale_y_continuous(name = paste("Média móvel (10 valores) de Beta para t=", as.character(i))) +
theme_bw()
}
par(old.par)
##############
# Volatilidade
##############
# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
plot(mcmc_sd_selic[,i], auto.layout = T, main = paste("D.P. em t=", as.character(i), "da Selic"))
}
par(old.par)
# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
df1 <- data.frame(sequencia2,sqrt(mcmc_sd_selic[,i]))
names(df1) <- c("Iteracao", "DP")
ggplot(df1, aes(Iteracao, DP)) +
geom_line(aes(y = rollmean(DP, 10, na.pad = TRUE))) +
scale_y_continuous(name = paste("Média móvel (10 valores) do desvio padrão para t=", as.character(i))) +
theme_bw()
}
par(old.par)
#####################################
# Para o IBC-Br
#####################################
#############
# Betas
#############
# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
plot(mcmc_beta_ibcbr[,i], auto.layout = T, main = paste("Beta em t=", as.character(i), "do IBC-Br"))
}
par(old.par)
# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
df1 <- data.frame(sequencia2,mcmc_beta_ibcbr[,i])
names(df1) <- c("Iteracao", "Beta")
ggplot(df1, aes(Iteracao, Beta)) +
geom_line(aes(y = rollmean(Beta, 10, na.pad = TRUE))) +
scale_y_continuous(name = paste("Média móvel (10 valores) de Beta para t=", as.character(i))) +
theme_bw()
}
par(old.par)
##############
# Volatilidade
##############
# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
plot(mcmc_sd_ibcbr[,i], auto.layout = T, main = paste("D.P. em t=", as.character(i), "da Selic"))
}
par(old.par)
# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
df1 <- data.frame(sequencia2,sqrt(mcmc_sd_ibcbr[,i]))
names(df1) <- c("Iteracao", "DP")
ggplot(df1, aes(Iteracao, DP)) +
geom_line(aes(y = rollmean(DP, 10, na.pad = TRUE))) +
scale_y_continuous(name = paste("Média móvel (10 valores) do desvio padrão para t=", as.character(i))) +
theme_bw()
}
par(old.par)
#####################################
# Para o Câmbio
#####################################
#############
# Betas
#############
# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
plot(mcmc_beta_cambio[,i], auto.layout = T, main = paste("Beta em t=", as.character(i), "do Câmbio"))
}
par(old.par)
# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
df1 <- data.frame(sequencia2,mcmc_beta_cambio[,i])
names(df1) <- c("Iteracao", "Beta")
ggplot(df1, aes(Iteracao, Beta)) +
geom_line(aes(y = rollmean(Beta, 10, na.pad = TRUE))) +
scale_y_continuous(name = paste("Média móvel (10 valores) de Beta para t=", as.character(i))) +
theme_bw()
}
par(old.par)
##############
# Volatilidade
##############
# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
plot(mcmc_sd_cambio[,i], auto.layout = T, main = paste("D.P. em t=", as.character(i), "da Selic"))
}
par(old.par)
# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
df1 <- data.frame(sequencia2,sqrt(mcmc_sd_cambio[,i]))
names(df1) <- c("Iteracao", "DP")
ggplot(df1, aes(Iteracao, DP)) +
geom_line(aes(y = rollmean(DP, 10, na.pad = TRUE))) +
scale_y_continuous(name = paste("Média móvel (10 valores) do desvio padrão para t=", as.character(i))) +
theme_bw()
}
par(old.par)
#####################################
# Para a IPCA
#####################################
#############
# Betas
#############
# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
plot(mcmc_beta_ipca[,i], auto.layout = T, main = paste("Beta em t=", as.character(i), "do IPCA"))
}
par(old.par)
# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
df1 <- data.frame(sequencia2,mcmc_beta_ipca[,i])
names(df1) <- c("Iteracao", "Beta")
ggplot(df1, aes(Iteracao, Beta)) +
geom_line(aes(y = rollmean(Beta, 10, na.pad = TRUE))) +
scale_y_continuous(name = paste("Média móvel (10 valores) de Beta para t=", as.character(i))) +
theme_bw()
}
par(old.par)
##############
# Volatilidade
##############
# Traceplot e density plots
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
plot(mcmc_sd_ipca[,i], auto.layout = T, main = paste("D.P. em t=", as.character(i), "do IPCA"))
}
par(old.par)
# Gráfico da média com janela móvel
old.par <- par(mfrow=c(5,3))
for (i in sequencia){
df1 <- data.frame(sequencia2,sqrt(mcmc_sd_ipca[,i]))
names(df1) <- c("Iteracao", "DP")
ggplot(df1, aes(Iteracao, DP)) +
geom_line(aes(y = rollmean(DP, 10, na.pad = TRUE))) +
scale_y_continuous(name = paste("Média móvel (10 valores) do desvio padrão para t=", as.character(i))) +
theme_bw()
}
par(old.par)
Agora tem o cumuplot() que faz o gráfico dos quantis da cadeia (de forma acumulada). O default é 0.025, 0.5, 0.975. Se eu tiver tempo depois, vou mexer no layout desses gráficos pois são feiosos demais.
#####################################
# Para a razão capital trabalho
#####################################
#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
cumuplot(mcmc_beta_ct[1:50000,i], auto.layout = F, main = paste("Quantis de Beta em t=", as.character(i), "da Razão Cap./Trab."))
}
#par(old.par)
#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
cumuplot(mcmc_sd_ct[1:50000,i], auto.layout = F, main = paste("Quantis do D.P. em t=", as.character(i), "da Razão Cap./Trab."))
}
#par(old.par)
#####################################
# Para a selic
#####################################
#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
cumuplot(mcmc_beta_selic[1:50000,i], auto.layout = F, main = paste("Quantis de Beta em t=", as.character(i), "da Selic"))
}
#par(old.par)
#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
cumuplot(mcmc_sd_selic[1:50000,i], auto.layout = F, main = paste("Quantis do D.P. em t=", as.character(i), "da Selic"))
}
#par(old.par)
#####################################
# Para o IBC-Br
#####################################
#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
cumuplot(mcmc_beta_ibcbr[1:50000,i], auto.layout = F, main = paste("Quantis de Beta em t=", as.character(i), "do IBC-Br"))
}
#par(old.par)
#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
cumuplot(mcmc_sd_ibcbr[1:50000,i], auto.layout = F, main = paste("Quantis do D.P. em t=", as.character(i), "do IBC-Br"))
}
#par(old.par)
#####################################
# Para o Câmbio
#####################################
#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
cumuplot(mcmc_beta_cambio[1:50000,i], auto.layout = F, main = paste("Quantis de Beta em t=", as.character(i), "do Câmbio"))
}
#par(old.par)
#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
cumuplot(mcmc_sd_cambio[1:50000,i], auto.layout = F, main = paste("Quantis do D.P. em t=", as.character(i), "do Câmbio"))
}
#par(old.par)
#####################################
# Para o IPCA
#####################################
#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
cumuplot(mcmc_beta_ipca[1:50000,i], auto.layout = F, main = paste("Quantis de Beta em t=", as.character(i), "do IPCA"))
}
#par(old.par)
#old.par <- par(mfrow=c(5,3))
for (i in sequencia){
cumuplot(mcmc_sd_ipca[1:50000,i], auto.layout = F, main = paste("Quantis do D.P. em t=", as.character(i), "do IPCA"))
}
#par(old.par)
Por fim, tem o diagnóstico do Geweke, que calcula uma estatística com distribuição normal e que compara o início e o final da distribuição. De novo vou plotar só uma parte para não ficar muito grande
geweke.plot(mcmc_beta_ct[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)
geweke.plot(mcmc_beta_selic[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)
geweke.plot(mcmc_beta_ibcbr[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)
geweke.plot(mcmc_beta_cambio[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)
geweke.plot(mcmc_beta_ipca[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)
geweke.plot(mcmc_sd_ct[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)
geweke.plot(mcmc_sd_selic[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)
geweke.plot(mcmc_sd_ibcbr[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)
geweke.plot(mcmc_sd_cambio[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)
geweke.plot(mcmc_sd_ipca[,sequencia], frac1 = .10, frac2 = .10, pvalue = 0.05)
O bloco abaixo eu montei para entender melhor o que a autocorr() do pacote coda faz. No fim, vou ficar com a função de autocorrelação (não a parcial) pois ela avalia o efeito até o instante j.
# Tentando descobrir o que as coisas são
# Esse é o do coda, ele é o mais parecido com a acf, mas não chega a ser igual.
autocorr(mcmc_sd_ct[,i], lags = c(20,100,1000))
c(acf(mcmc_sd_ct[,1], lag.max = 20, plot = F)[20]$acf,acf(mcmc_sd_ct[,1], lag.max = 100, plot = F)[100]$acf, acf(mcmc_sd_ct[,1], lag.max = 1000, plot = F)[1000]$acf)
c(Acf(mcmc_sd_ct[,1], lag.max = 20, plot = F, type = "correlation")[20]$acf, Acf(mcmc_sd_ct[,1], lag.max = 100, plot = F, type = "correlation")[100]$acf, Acf(mcmc_sd_ct[,1], lag.max = 1000, plot = F, type = "correlation")[1000]$acf)
# Esse não é parecido com ninguém
c(Acf(mcmc_sd_ct[,1], lag.max = 20, plot = F, type = "covariance")[20]$acf, Acf(mcmc_sd_ct[,1], lag.max = 100, plot = F, type = "covariance")[100]$acf, Acf(mcmc_sd_ct[,1], lag.max = 1000, plot = F, type = "covariance")[1000]$acf)
# Esses 3 são iguais
# Function Pacf computes (and by default plots) an estimate of the partial autocorrelation function of a (possibly multivariate) time series
# Aviso sobre as funções Acf e Pacf: The functions improve the acf, pacf and ccf functions. The main differences are that Acf does not plot a spike at lag 0 when type=="correlation" (which is redundant) and the horizontal axes show lags in time units rather than seasonal units..
c(pacf(mcmc_sd_ct[,1], lag.max = 20, plot = F)[20]$acf,pacf(mcmc_sd_ct[,1], lag.max = 100, plot = F)[100]$acf, pacf(mcmc_sd_ct[,1], lag.max = 1000, plot = F)[1000]$acf)
c(Acf(mcmc_sd_ct[,1], lag.max = 20, plot = F, type = "partial")[20]$acf, Acf(mcmc_sd_ct[,1], lag.max = 100, plot = F, type = "partial")[100]$acf, Acf(mcmc_sd_ct[,1], lag.max = 1000, plot = F, type = "partial")[1000]$acf)
c(Pacf(mcmc_sd_ct[,1], lag.max = 20, plot = F)[20]$acf, Pacf(mcmc_sd_ct[,1], lag.max = 100, plot = F)[100]$acf, Pacf(mcmc_sd_ct[,1], lag.max = 1000, plot = F)[1000]$acf)
Agora estou calculando a função te autocorrelação para depois fazer uns fancy graphs.
autocorrelacao_beta_ct <- vector()
autocorrelacao_beta_selic <- vector()
autocorrelacao_beta_ibcbr <- vector()
autocorrelacao_beta_cambio <- vector()
autocorrelacao_beta_ipca <- vector()
autocorrelacao_sd_ct <- vector()
autocorrelacao_sd_selic <- vector()
autocorrelacao_sd_ibcbr <- vector()
autocorrelacao_sd_cambio <- vector()
autocorrelacao_sd_ipca <- vector()
pm <- proc.time()
for (i in 1:dim(mcmc_beta_ct)[2]) {
autocorrelacao_beta_ct <- cbind(autocorrelacao_beta_ct, autocorr(mcmc_beta_ct[,i], lags = c(20,100,1000)))
autocorrelacao_beta_selic <- cbind(autocorrelacao_beta_selic, autocorr(mcmc_beta_selic[,i], lags = c(20,100,1000)))
autocorrelacao_beta_ibcbr <- cbind(autocorrelacao_beta_ibcbr, autocorr(mcmc_beta_ibcbr[,i], lags = c(20,100,1000)))
autocorrelacao_beta_cambio <- cbind(autocorrelacao_beta_cambio, autocorr(mcmc_beta_cambio[,i], lags = c(20,100,1000)))
autocorrelacao_beta_ipca <- cbind(autocorrelacao_beta_ipca, autocorr(mcmc_sd_ipca[,i], lags = c(20,100,1000)))
autocorrelacao_sd_ct <- cbind(autocorrelacao_sd_ct, autocorr(mcmc_sd_ct[,i], lags = c(20,100,1000)))
autocorrelacao_sd_selic <- cbind(autocorrelacao_sd_selic, autocorr(mcmc_sd_selic[,i], lags = c(20,100,1000)))
autocorrelacao_sd_ibcbr <- cbind(autocorrelacao_sd_ibcbr, autocorr(mcmc_sd_ibcbr[,i], lags = c(20,100,1000)))
autocorrelacao_sd_cambio <- cbind(autocorrelacao_sd_cambio, autocorr(mcmc_sd_cambio[,i], lags = c(20,100,1000)))
autocorrelacao_sd_ipca <- cbind(autocorrelacao_sd_ipca, autocorr(mcmc_sd_ipca[,i], lags = c(20,100,1000)))
}
proc.time() - pm
## user system elapsed
## 150.36 0.25 152.46
Agora vamos para os gráficos…. a ideia é fazer algo como o Primiceri tem no anexo B (figura 9 - topo), porém dividindo por variável do VAR.
# Montando os data frames
#############################################
#############################################
# Betas
#############################################
#############################################
autocorrelacao_beta_ct <- data.frame(1:152, t(autocorrelacao_beta_ct))
names(autocorrelacao_beta_ct) <- c("tempo", "20", "100", "1000")
autocorrelacao_beta_ct <- melt(autocorrelacao_beta_ct, id.vars = "tempo")
autocorrelacao_beta_ct[,2] <- factor(autocorrelacao_beta_ct[,2])
autocorrelacao_beta_selic1 <- data.frame(1:152, t(autocorrelacao_beta_selic))
names(autocorrelacao_beta_selic1) <- c("tempo", "20", "100", "1000")
autocorrelacao_beta_selic1 <- melt(autocorrelacao_beta_selic1, id.vars = "tempo")
autocorrelacao_beta_selic1[,2] <- factor(autocorrelacao_beta_selic1[,2])
autocorrelacao_beta_ibcbr <- data.frame(1:152, t(autocorrelacao_beta_ibcbr))
names(autocorrelacao_beta_ibcbr) <- c("tempo", "20", "100", "1000")
autocorrelacao_beta_ibcbr <- melt(autocorrelacao_beta_ibcbr, id.vars = "tempo")
autocorrelacao_beta_ibcbr[,2] <- factor(autocorrelacao_beta_ibcbr[,2])
autocorrelacao_beta_cambio <- data.frame(1:152, t(autocorrelacao_beta_cambio))
names(autocorrelacao_beta_cambio) <- c("tempo", "20", "100", "1000")
autocorrelacao_beta_cambio <- melt(autocorrelacao_beta_cambio, id.vars = "tempo")
autocorrelacao_beta_cambio[,2] <- factor(autocorrelacao_beta_cambio[,2])
autocorrelacao_beta_ipca <- data.frame(1:152, t(autocorrelacao_beta_ipca))
names(autocorrelacao_beta_ipca) <- c("tempo", "20", "100", "1000")
autocorrelacao_beta_ipca <- melt(autocorrelacao_beta_ipca, id.vars = "tempo")
autocorrelacao_beta_ipca[,2] <- factor(autocorrelacao_beta_ipca[,2])
#############################################
#############################################
# Volatilidade
#############################################
#############################################
autocorrelacao_sd_ct <- data.frame(1:152, t(autocorrelacao_sd_ct))
names(autocorrelacao_sd_ct) <- c("tempo", "20", "100", "1000")
autocorrelacao_sd_ct <- melt(autocorrelacao_sd_ct, id.vars = "tempo")
autocorrelacao_sd_ct[,2] <- factor(autocorrelacao_sd_ct[,2])
autocorrelacao_sd_selic <- data.frame(1:152, t(autocorrelacao_sd_selic))
names(autocorrelacao_sd_selic) <- c("tempo", "20", "100", "1000")
autocorrelacao_sd_selic <- melt(autocorrelacao_sd_selic, id.vars = "tempo")
autocorrelacao_sd_selic[,2] <- factor(autocorrelacao_sd_selic[,2])
autocorrelacao_sd_ibcbr <- data.frame(1:152, t(autocorrelacao_sd_ibcbr))
names(autocorrelacao_sd_ibcbr) <- c("tempo", "20", "100", "1000")
autocorrelacao_sd_ibcbr <- melt(autocorrelacao_sd_ibcbr, id.vars = "tempo")
autocorrelacao_sd_ibcbr[,2] <- factor(autocorrelacao_sd_ibcbr[,2])
autocorrelacao_sd_cambio <- data.frame(1:152, t(autocorrelacao_sd_cambio))
names(autocorrelacao_sd_cambio) <- c("tempo", "20", "100", "1000")
autocorrelacao_sd_cambio <- melt(autocorrelacao_sd_cambio, id.vars = "tempo")
autocorrelacao_sd_cambio[,2] <- factor(autocorrelacao_sd_cambio[,2])
autocorrelacao_sd_ipca <- data.frame(1:152, t(autocorrelacao_sd_ipca))
names(autocorrelacao_sd_ipca) <- c("tempo", "20", "100", "1000")
autocorrelacao_sd_ipca <- melt(autocorrelacao_sd_ipca, id.vars = "tempo")
autocorrelacao_sd_ipca[,2] <- factor(autocorrelacao_sd_ipca[,2])
#############################################
#############################################
# Gráficos
#############################################
#############################################
##################### BETAS ########################
########################
# Razão capital trabalho
########################
p1a <- ggplot(autocorrelacao_beta_ct, aes(tempo, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
scale_y_continuous(name="Razão Cap./Trab.", lim = c(-0.1,0.6)) +
scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
theme_bw() +
scale_colour_brewer(palette="Set1")
p1a <- p1a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", legend.text = element_text(size = 8), legend.title = element_text(size = 8), axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.05,0.1,0.1), "cm"))
p1a <- p1a + labs(linetype = "Lags") # I couldn't fix the legend name properly, so why not use a gambiarra?
p1a <- p1a + labs(colour='Lags')
########################
# selic
########################
p2a <- ggplot(autocorrelacao_beta_selic1, aes(tempo, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
scale_y_continuous(name="Selic", lim = c(-0.1,0.6)) +
scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
theme_bw() +
scale_colour_brewer(palette="Set1")
p2a <- p2a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.05,0.1,0.1), "cm"))
p2a <- p2a + labs(linetype = "Lags")
p2a <- p2a + labs(colour='Lags')
########################
# ibcbr
########################
p3a <- ggplot(autocorrelacao_beta_ibcbr, aes(tempo, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
scale_y_continuous(name="IBCBr", lim = c(-0.1,0.6)) +
scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
theme_bw() +
scale_colour_brewer(palette="Set1")
p3a <- p3a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="right", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.05,0.1,0.1), "cm"))
p3a <- p3a + labs(linetype = "Lags")
p3a <- p3a + labs(colour='Lags')
########################
# cambio
########################
p4a <- ggplot(autocorrelacao_beta_cambio, aes(tempo, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
scale_y_continuous(name="Câmbio", lim = c(-0.1,0.6)) +
scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
theme_bw() +
scale_colour_brewer(palette="Set1")
p4a <- p4a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.05,0.1,0.1), "cm"))
p4a <- p4a + labs(linetype = "Lags")
p4a <- p4a + labs(colour='Lags')
########################
# ipca
########################
p5a <- ggplot(autocorrelacao_beta_ipca, aes(tempo, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
scale_y_continuous(name="IPCA", lim = c(-0.1,0.6)) +
scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
theme_bw() +
scale_colour_brewer(palette="Set1")
p5a <- p5a + theme(axis.text.x = element_text(angle=25, vjust = 2, size = 5, hjust = 1), axis.title.x = element_text(size = 6, vjust = 5), axis.title.y = element_text(size =6), legend.position="right", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.05,0.0001,0.1), "cm")) # cima,direita, baixo, esquerda)
p5a <- p5a + labs(linetype = "Lags") # I couldn't fix the legend name properly, so why not?
p5a <- p5a + labs(colour='Lags')
#grid.draw(rbind(ggplotGrob(p1a), ggplotGrob(p2a), ggplotGrob(p3a), ggplotGrob(p4a), ggplotGrob(p5a), size = "first"))
################## VOLATILIDADE ######################
########################
# Razão capital trabalho
########################
p1 <- ggplot(autocorrelacao_sd_ct, aes(tempo, value, colour = variable)) +
geom_line(alpha = 1, show.legend=T, aes(linetype = variable))+
scale_y_continuous(name="Razão Cap./Trab.", lim = c(-0.1,0.5)) +
scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
theme_bw() +
scale_colour_brewer(palette="Set1")
p1 <- p1 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", legend.text = element_text(size = 8), legend.title = element_text(size = 8), axis.text.y = element_text(size = 6), legend.key.size = unit(.5, "cm"), plot.margin = unit(c(0.009,0.08,0.1,0.1), "cm"))
p1 <- p1 + labs(linetype = "Lags") # I couldn't fix the legend name properly, so why not use a gambiarra?
p1 <- p1 + labs(colour='Lags')
########################
# selic
########################
p2 <- ggplot(autocorrelacao_sd_selic, aes(tempo, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
scale_y_continuous(name="Selic", lim = c(-0.1,0.5)) +
scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
theme_bw() +
scale_colour_brewer(palette="Set1")
p2 <- p2 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.08,0.1,0.1), "cm"))
p2 <- p2 + labs(linetype = "Lags")
p2 <- p2 + labs(colour='Lags')
########################
# ibcbr
########################
p3 <- ggplot(autocorrelacao_sd_ibcbr, aes(tempo, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
scale_y_continuous(name="IBCBr", lim = c(-0.1,0.5)) +
scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
theme_bw() +
scale_colour_brewer(palette="Set1")
p3 <- p3 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="right", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.08,0.1,0.1), "cm"))
p3 <- p3 + labs(linetype = "Lags")
p3 <- p3 + labs(colour='Lags')
########################
# cambio
########################
p4 <- ggplot(autocorrelacao_sd_cambio, aes(tempo, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
scale_y_continuous(name="Câmbio", lim = c(-0.1,0.5)) +
scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
theme_bw() +
scale_colour_brewer(palette="Set1")
p4 <- p4 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6), legend.position="top", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.08,0.1,0.1), "cm")) # cima,direita, baixo, esquerda)
p4 <- p4 + labs(linetype = "Lags")
p4 <- p4 + labs(colour='Lags')
########################
# ipca
########################
p5 <- ggplot(autocorrelacao_sd_ipca, aes(tempo, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, aes(linetype = variable))+
scale_y_continuous(name="IPCA", lim = c(-0.1,0.5)) +
scale_x_continuous(name = "Tempo", breaks = seq(0,155,10), expand = c(0,0)) +
theme_bw() +
scale_colour_brewer(palette="Set1")
p5 <- p5 + theme(axis.text.x = element_text(angle=25, vjust = 2, size = 5, hjust = 1), axis.title.x = element_text(size = 6, vjust = 5), axis.title.y = element_text(size =6), legend.position="right", axis.text.y = element_text(size = 6), plot.margin = unit(c(0.009,0.08,0.0001,0.1), "cm")) # cima,direita, baixo, esquerda
p5 <- p5 + labs(linetype = "Lags") # I couldn't fix the legend name properly, so why not?
p5 <- p5 + labs(colour='Lags')
ggarrange(p1a, p1, p2a, p2, p3a, p3, p4a, p4, p5a, p5, common.legend = TRUE, ncol = 2, nrow = 5)
#################
# grid.draw(rbind(ggplotGrob(p1), ggplotGrob(p2), ggplotGrob(p3), ggplotGrob(p4), ggplotGrob(p5), size = "first"))
# grid.arrange(p1, p2, p3, p4, p5, ncol=1, nrow = 5) # This squish the graph with the legend
# Estimar outroS 4 modelos guardando apenas as últimas 200 observações e fazer uma anova não paramétrica
A estatística de Gelman-Rubin mede se há diferença significativa entre a variância dentro e entre diversas cadeias utilizando um fator de redução de escala. Aqui eu optei por simular duas cadeias e fazer a comparação entre elas. Em função da memória do computador eu utilizei um número menor de iterações salvas: ao invés de salvar todos os valores, armazenei apenas os 10.000 últimos valores (20%).
# Primeira coisa é limpar a memória
rm(list = ls())
# Auxiliary variables, so I don't need to bother when something changes
inicio <- "2003-02-01"
fim <- "2017-10-31"
# Don't mess with this code
inicio_cambio <- paste(seq(as.Date(inicio), length = 2, by = "-1 month")[2]) # 1 mês antes do início das outras séries
inicio_ipca <- paste(seq(as.Date(inicio), length = 12, by = "-1 month")[12]) # 12 meses antes do início das outras séries
trabalho <- BETS.get("7620", from = inicio, to = fim)
capital <- BETS.get("7621", from = inicio, to = fim)
capital_trabalho <- capital/trabalho
# Usando ndiffs(ibcbr,test="adf",alpha = 0.1) se encontra que o ibcbr é não estacionário.
# Como todas as outras séries são estacionárias, eu vou fazer também a diferença do log, igual no câmbio
ibcbr_raw <- BETS.get("24364", from = inicio_cambio, to = fim)
ibcbr <- diff(log(ibcbr_raw), 1)
cambio_raw <- BETS.get("3696", from = inicio_cambio, to = fim)
cambio <- diff(log(cambio_raw), 1)
selic_4390 <- BETS.get("4390", from = inicio, to = fim)
selic <- ((1+selic_4390/100)^(12)-1)*100
ipca_raw <- BETS.get("433", from = inicio_ipca, to = fim)
ipca_acum <- ipca_raw/100 + 1
ipca <- vector()
final <- length(ipca_acum)
for (i in 12:final){
ipca[(i-11)] <- (prod(ipca_acum[(i-11):i])-1)*100
}
ano <- as.numeric(substr(inicio, start = 1, stop = 4))
mes <- as.numeric(substr(inicio, start = 6, stop = 7))
dia <- as.numeric(substr(inicio, start = 9, stop = 10))
ipca <- ts(ipca, start = c(ano, mes, dia), frequency = 12) # Date format YYYY MM DD
m <- seas(x = capital_trabalho)
capital_trabalho2 <- final(m)
m <- seas(x = selic, transform.function = "none")
selic2 <- final(m)
var1 <- cbind(capital_trabalho2, selic2, ibcbr, cambio, ipca) # The bvar function does not allows data.frames
nburn. <- 5000
nrep. <- 50000
thinfac. <- 5
set.seed(6969)
pm <- proc.time()
fit1 <- bvar.sv.tvp(var1, p=1, tau = 24, nburn = nburn., nrep = nrep., nf = 1, thinfac = thinfac.)
## [1] "2018-01-26 10:01:39 -- now starting MCMC"
## [1] "2018-01-26 10:02:48 -- now at iteration 5000"
## [1] "2018-01-26 10:05:04 -- now at iteration 15000"
## [1] "2018-01-26 10:07:20 -- now at iteration 25000"
## [1] "2018-01-26 10:09:36 -- now at iteration 35000"
## [1] "2018-01-26 10:11:53 -- now at iteration 45000"
## [1] "2018-01-26 10:14:09 -- now at iteration 55000"
proc.time() - pm
## user system elapsed
## 733.25 16.95 750.40
set.seed(9696)
fit2 <- bvar.sv.tvp(var1, p=1, tau = 24, nburn = nburn., nrep = nrep., nf = 1, thinfac = thinfac.)
## [1] "2018-01-26 10:14:10 -- now starting MCMC"
## [1] "2018-01-26 10:15:17 -- now at iteration 5000"
## [1] "2018-01-26 10:17:34 -- now at iteration 15000"
## [1] "2018-01-26 10:19:50 -- now at iteration 25000"
## [1] "2018-01-26 10:22:06 -- now at iteration 35000"
## [1] "2018-01-26 10:24:22 -- now at iteration 45000"
## [1] "2018-01-26 10:26:40 -- now at iteration 55000"
mcmc_1 <- mcmc(cbind(parameter.draws(fit1, type = "lag1", row = 1, col = 1), parameter.draws(fit1, type = "lag1", row = 2, col = 2), parameter.draws(fit1, type = "lag1", row = 3, col = 3), parameter.draws(fit1, type = "lag1", row = 4, col = 4), parameter.draws(fit1, type = "lag1", row = 5, col = 5)))
mcmc_2 <- mcmc(cbind(parameter.draws(fit2, type = "lag1", row = 1, col = 1), parameter.draws(fit2, type = "lag1", row = 2, col = 2), parameter.draws(fit2, type = "lag1", row = 3, col = 3), parameter.draws(fit2, type = "lag1", row = 4, col = 4), parameter.draws(fit2, type = "lag1", row = 5, col = 5)))
lista <- mcmc.list(mcmc_1, mcmc_2)
gelman.diag(lista)
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## [1,] 1.01 1.03
## [2,] 1.01 1.03
## [3,] 1.01 1.03
## [4,] 1.01 1.03
## [5,] 1.01 1.02
## [6,] 1.01 1.03
## [7,] 1.00 1.02
## [8,] 1.00 1.02
## [9,] 1.00 1.02
## [10,] 1.00 1.02
## [11,] 1.00 1.02
## [12,] 1.00 1.01
## [13,] 1.00 1.01
## [14,] 1.00 1.01
## [15,] 1.00 1.01
## [16,] 1.00 1.01
## [17,] 1.00 1.02
## [18,] 1.00 1.01
## [19,] 1.00 1.01
## [20,] 1.00 1.01
## [21,] 1.00 1.01
## [22,] 1.00 1.01
## [23,] 1.00 1.01
## [24,] 1.00 1.02
## [25,] 1.01 1.02
## [26,] 1.01 1.03
## [27,] 1.01 1.03
## [28,] 1.01 1.03
## [29,] 1.01 1.03
## [30,] 1.00 1.02
## [31,] 1.00 1.01
## [32,] 1.00 1.01
## [33,] 1.00 1.01
## [34,] 1.00 1.01
## [35,] 1.00 1.01
## [36,] 1.00 1.01
## [37,] 1.00 1.01
## [38,] 1.00 1.02
## [39,] 1.01 1.03
## [40,] 1.01 1.03
## [41,] 1.00 1.02
## [42,] 1.00 1.01
## [43,] 1.00 1.01
## [44,] 1.00 1.01
## [45,] 1.00 1.01
## [46,] 1.00 1.00
## [47,] 1.00 1.00
## [48,] 1.00 1.00
## [49,] 1.00 1.00
## [50,] 1.00 1.00
## [51,] 1.00 1.01
## [52,] 1.00 1.01
## [53,] 1.00 1.01
## [54,] 1.00 1.00
## [55,] 1.00 1.00
## [56,] 1.00 1.00
## [57,] 1.00 1.00
## [58,] 1.00 1.00
## [59,] 1.00 1.00
## [60,] 1.00 1.01
## [61,] 1.00 1.01
## [62,] 1.00 1.00
## [63,] 1.00 1.00
## [64,] 1.00 1.00
## [65,] 1.00 1.00
## [66,] 1.00 1.01
## [67,] 1.00 1.02
## [68,] 1.01 1.04
## [69,] 1.01 1.05
## [70,] 1.01 1.04
## [71,] 1.01 1.04
## [72,] 1.01 1.04
## [73,] 1.01 1.04
## [74,] 1.01 1.04
## [75,] 1.01 1.04
## [76,] 1.01 1.04
## [77,] 1.01 1.04
## [78,] 1.01 1.03
## [79,] 1.00 1.02
## [80,] 1.00 1.00
## [81,] 1.00 1.00
## [82,] 1.00 1.00
## [83,] 1.00 1.00
## [84,] 1.00 1.00
## [85,] 1.00 1.00
## [86,] 1.00 1.00
## [87,] 1.00 1.00
## [88,] 1.00 1.00
## [89,] 1.00 1.01
## [90,] 1.00 1.00
## [91,] 1.00 1.00
## [92,] 1.00 1.00
## [93,] 1.00 1.00
## [94,] 1.00 1.00
## [95,] 1.00 1.00
## [96,] 1.00 1.00
## [97,] 1.00 1.00
## [98,] 1.00 1.00
## [99,] 1.00 1.00
## [100,] 1.00 1.00
## [101,] 1.00 1.01
## [102,] 1.00 1.01
## [103,] 1.00 1.01
## [104,] 1.00 1.00
## [105,] 1.00 1.00
## [106,] 1.00 1.00
## [107,] 1.00 1.00
## [108,] 1.00 1.00
## [109,] 1.00 1.00
## [110,] 1.00 1.00
## [111,] 1.00 1.01
## [112,] 1.00 1.01
## [113,] 1.00 1.01
## [114,] 1.00 1.01
## [115,] 1.00 1.01
## [116,] 1.00 1.01
## [117,] 1.00 1.02
## [118,] 1.01 1.03
## [119,] 1.01 1.06
## [120,] 1.01 1.07
## [121,] 1.02 1.07
## [122,] 1.02 1.07
## [123,] 1.02 1.07
## [124,] 1.02 1.09
## [125,] 1.02 1.09
## [126,] 1.02 1.08
## [127,] 1.02 1.07
## [128,] 1.02 1.07
## [129,] 1.01 1.06
## [130,] 1.01 1.03
## [131,] 1.00 1.01
## [132,] 1.00 1.00
## [133,] 1.00 1.00
## [134,] 1.00 1.00
## [135,] 1.00 1.01
## [136,] 1.00 1.02
## [137,] 1.01 1.03
## [138,] 1.01 1.04
## [139,] 1.02 1.08
## [140,] 1.03 1.11
## [141,] 1.03 1.13
## [142,] 1.03 1.15
## [143,] 1.04 1.15
## [144,] 1.03 1.15
## [145,] 1.03 1.14
## [146,] 1.03 1.14
## [147,] 1.03 1.14
## [148,] 1.03 1.13
## [149,] 1.03 1.12
## [150,] 1.02 1.10
## [151,] 1.02 1.08
## [152,] 1.02 1.06
## [153,] 1.01 1.03
## [154,] 1.01 1.03
## [155,] 1.00 1.01
## [156,] 1.00 1.00
## [157,] 1.00 1.01
## [158,] 1.00 1.00
## [159,] 1.00 1.00
## [160,] 1.00 1.01
## [161,] 1.00 1.00
## [162,] 1.00 1.00
## [163,] 1.00 1.01
## [164,] 1.01 1.03
## [165,] 1.01 1.04
## [166,] 1.03 1.10
## [167,] 1.03 1.09
## [168,] 1.03 1.09
## [169,] 1.02 1.07
## [170,] 1.02 1.06
## [171,] 1.02 1.06
## [172,] 1.02 1.06
## [173,] 1.02 1.05
## [174,] 1.02 1.05
## [175,] 1.02 1.05
## [176,] 1.02 1.05
## [177,] 1.02 1.05
## [178,] 1.02 1.05
## [179,] 1.01 1.05
## [180,] 1.02 1.05
## [181,] 1.02 1.06
## [182,] 1.02 1.05
## [183,] 1.02 1.05
## [184,] 1.02 1.05
## [185,] 1.02 1.05
## [186,] 1.02 1.06
## [187,] 1.02 1.07
## [188,] 1.02 1.07
## [189,] 1.02 1.07
## [190,] 1.02 1.09
## [191,] 1.02 1.09
## [192,] 1.02 1.09
## [193,] 1.02 1.06
## [194,] 1.01 1.05
## [195,] 1.01 1.04
## [196,] 1.01 1.02
## [197,] 1.01 1.01
## [198,] 1.00 1.00
## [199,] 1.00 1.00
## [200,] 1.00 1.00
## [201,] 1.01 1.01
## [202,] 1.01 1.02
## [203,] 1.01 1.02
## [204,] 1.01 1.03
## [205,] 1.01 1.02
## [206,] 1.01 1.01
## [207,] 1.01 1.01
## [208,] 1.00 1.00
## [209,] 1.00 1.01
## [210,] 1.01 1.01
## [211,] 1.01 1.02
## [212,] 1.01 1.03
## [213,] 1.01 1.03
## [214,] 1.01 1.04
## [215,] 1.01 1.04
## [216,] 1.01 1.03
## [217,] 1.01 1.04
## [218,] 1.01 1.05
## [219,] 1.02 1.07
## [220,] 1.02 1.09
## [221,] 1.02 1.09
## [222,] 1.03 1.10
## [223,] 1.03 1.11
## [224,] 1.03 1.10
## [225,] 1.03 1.10
## [226,] 1.02 1.09
## [227,] 1.02 1.08
## [228,] 1.02 1.06
## [229,] 1.02 1.05
## [230,] 1.01 1.03
## [231,] 1.00 1.01
## [232,] 1.00 1.00
## [233,] 1.00 1.01
## [234,] 1.01 1.01
## [235,] 1.01 1.02
## [236,] 1.01 1.03
## [237,] 1.01 1.03
## [238,] 1.01 1.03
## [239,] 1.01 1.03
## [240,] 1.01 1.02
## [241,] 1.00 1.00
## [242,] 1.00 1.00
## [243,] 1.00 1.00
## [244,] 1.00 1.00
## [245,] 1.00 1.01
## [246,] 1.01 1.02
## [247,] 1.01 1.02
## [248,] 1.01 1.03
## [249,] 1.01 1.02
## [250,] 1.01 1.02
## [251,] 1.01 1.02
## [252,] 1.02 1.08
## [253,] 1.01 1.02
## [254,] 1.01 1.03
## [255,] 1.01 1.02
## [256,] 1.01 1.02
## [257,] 1.01 1.03
## [258,] 1.01 1.03
## [259,] 1.01 1.02
## [260,] 1.01 1.03
## [261,] 1.01 1.02
## [262,] 1.01 1.02
## [263,] 1.00 1.02
## [264,] 1.00 1.02
## [265,] 1.01 1.03
## [266,] 1.01 1.03
## [267,] 1.01 1.03
## [268,] 1.01 1.02
## [269,] 1.01 1.03
## [270,] 1.01 1.04
## [271,] 1.02 1.08
## [272,] 1.02 1.08
## [273,] 1.02 1.07
## [274,] 1.02 1.06
## [275,] 1.02 1.05
## [276,] 1.02 1.06
## [277,] 1.02 1.06
## [278,] 1.02 1.04
## [279,] 1.02 1.03
## [280,] 1.02 1.03
## [281,] 1.02 1.02
## [282,] 1.01 1.01
## [283,] 1.01 1.01
## [284,] 1.01 1.01
## [285,] 1.01 1.01
## [286,] 1.01 1.01
## [287,] 1.01 1.01
## [288,] 1.01 1.01
## [289,] 1.01 1.02
## [290,] 1.01 1.02
## [291,] 1.02 1.04
## [292,] 1.03 1.06
## [293,] 1.04 1.08
## [294,] 1.04 1.08
## [295,] 1.03 1.08
## [296,] 1.03 1.07
## [297,] 1.03 1.06
## [298,] 1.03 1.07
## [299,] 1.03 1.06
## [300,] 1.03 1.05
## [301,] 1.02 1.04
## [302,] 1.02 1.03
## [303,] 1.01 1.02
## [304,] 1.01 1.02
## [305,] 1.00 1.01
## [306,] 1.00 1.01
## [307,] 1.00 1.00
## [308,] 1.00 1.00
## [309,] 1.00 1.00
## [310,] 1.00 1.00
## [311,] 1.00 1.00
## [312,] 1.00 1.00
## [313,] 1.00 1.00
## [314,] 1.00 1.00
## [315,] 1.00 1.00
## [316,] 1.00 1.00
## [317,] 1.00 1.00
## [318,] 1.00 1.01
## [319,] 1.00 1.01
## [320,] 1.00 1.01
## [321,] 1.00 1.00
## [322,] 1.00 1.00
## [323,] 1.00 1.01
## [324,] 1.00 1.01
## [325,] 1.00 1.01
## [326,] 1.00 1.01
## [327,] 1.00 1.01
## [328,] 1.00 1.01
## [329,] 1.00 1.01
## [330,] 1.00 1.01
## [331,] 1.00 1.01
## [332,] 1.00 1.01
## [333,] 1.00 1.01
## [334,] 1.00 1.01
## [335,] 1.00 1.01
## [336,] 1.00 1.01
## [337,] 1.00 1.01
## [338,] 1.00 1.01
## [339,] 1.00 1.01
## [340,] 1.00 1.01
## [341,] 1.00 1.01
## [342,] 1.00 1.00
## [343,] 1.00 1.00
## [344,] 1.00 1.00
## [345,] 1.00 1.00
## [346,] 1.00 1.01
## [347,] 1.00 1.01
## [348,] 1.00 1.01
## [349,] 1.00 1.01
## [350,] 1.00 1.02
## [351,] 1.00 1.02
## [352,] 1.00 1.02
## [353,] 1.00 1.02
## [354,] 1.00 1.02
## [355,] 1.01 1.02
## [356,] 1.01 1.02
## [357,] 1.01 1.03
## [358,] 1.01 1.02
## [359,] 1.00 1.02
## [360,] 1.00 1.02
## [361,] 1.00 1.01
## [362,] 1.00 1.01
## [363,] 1.00 1.01
## [364,] 1.00 1.01
## [365,] 1.00 1.01
## [366,] 1.00 1.01
## [367,] 1.00 1.01
## [368,] 1.00 1.01
## [369,] 1.00 1.01
## [370,] 1.00 1.01
## [371,] 1.00 1.00
## [372,] 1.00 1.00
## [373,] 1.00 1.00
## [374,] 1.00 1.00
## [375,] 1.00 1.00
## [376,] 1.00 1.00
## [377,] 1.00 1.00
## [378,] 1.00 1.00
## [379,] 1.00 1.00
## [380,] 1.00 1.00
## [381,] 1.00 1.00
## [382,] 1.00 1.00
## [383,] 1.00 1.00
## [384,] 1.00 1.01
## [385,] 1.00 1.01
## [386,] 1.00 1.01
## [387,] 1.00 1.02
## [388,] 1.00 1.02
## [389,] 1.00 1.02
## [390,] 1.00 1.02
## [391,] 1.00 1.02
## [392,] 1.00 1.02
## [393,] 1.00 1.02
## [394,] 1.00 1.02
## [395,] 1.00 1.02
## [396,] 1.00 1.02
## [397,] 1.00 1.02
## [398,] 1.00 1.01
## [399,] 1.00 1.01
## [400,] 1.00 1.01
## [401,] 1.00 1.02
## [402,] 1.00 1.02
## [403,] 1.00 1.02
## [404,] 1.01 1.02
## [405,] 1.01 1.02
## [406,] 1.00 1.02
## [407,] 1.00 1.02
## [408,] 1.00 1.02
## [409,] 1.00 1.01
## [410,] 1.00 1.01
## [411,] 1.00 1.01
## [412,] 1.00 1.01
## [413,] 1.00 1.01
## [414,] 1.00 1.01
## [415,] 1.00 1.01
## [416,] 1.00 1.01
## [417,] 1.00 1.01
## [418,] 1.00 1.01
## [419,] 1.00 1.01
## [420,] 1.00 1.01
## [421,] 1.00 1.01
## [422,] 1.00 1.00
## [423,] 1.00 1.00
## [424,] 1.00 1.00
## [425,] 1.00 1.00
## [426,] 1.00 1.00
## [427,] 1.00 1.00
## [428,] 1.00 1.00
## [429,] 1.00 1.00
## [430,] 1.00 1.00
## [431,] 1.00 1.00
## [432,] 1.00 1.00
## [433,] 1.00 1.00
## [434,] 1.00 1.00
## [435,] 1.00 1.01
## [436,] 1.00 1.01
## [437,] 1.00 1.01
## [438,] 1.00 1.01
## [439,] 1.00 1.01
## [440,] 1.00 1.01
## [441,] 1.00 1.01
## [442,] 1.00 1.01
## [443,] 1.00 1.01
## [444,] 1.01 1.02
## [445,] 1.01 1.02
## [446,] 1.01 1.02
## [447,] 1.01 1.02
## [448,] 1.01 1.02
## [449,] 1.01 1.02
## [450,] 1.01 1.02
## [451,] 1.01 1.02
## [452,] 1.01 1.02
## [453,] 1.01 1.02
## [454,] 1.01 1.02
## [455,] 1.01 1.02
## [456,] 1.01 1.01
## [457,] 1.01 1.06
## [458,] 1.01 1.05
## [459,] 1.01 1.04
## [460,] 1.01 1.03
## [461,] 1.01 1.03
## [462,] 1.01 1.03
## [463,] 1.01 1.03
## [464,] 1.01 1.03
## [465,] 1.00 1.02
## [466,] 1.00 1.01
## [467,] 1.00 1.01
## [468,] 1.00 1.01
## [469,] 1.00 1.01
## [470,] 1.00 1.00
## [471,] 1.00 1.00
## [472,] 1.00 1.00
## [473,] 1.00 1.00
## [474,] 1.00 1.00
## [475,] 1.00 1.00
## [476,] 1.00 1.00
## [477,] 1.00 1.00
## [478,] 1.00 1.00
## [479,] 1.00 1.00
## [480,] 1.00 1.00
## [481,] 1.00 1.00
## [482,] 1.00 1.00
## [483,] 1.00 1.00
## [484,] 1.00 1.00
## [485,] 1.00 1.00
## [486,] 1.00 1.00
## [487,] 1.00 1.00
## [488,] 1.00 1.00
## [489,] 1.00 1.00
## [490,] 1.00 1.01
## [491,] 1.00 1.01
## [492,] 1.00 1.01
## [493,] 1.00 1.01
## [494,] 1.00 1.01
## [495,] 1.00 1.01
## [496,] 1.00 1.01
## [497,] 1.00 1.01
## [498,] 1.00 1.01
## [499,] 1.00 1.02
## [500,] 1.00 1.02
## [501,] 1.00 1.01
## [502,] 1.00 1.01
## [503,] 1.00 1.01
## [504,] 1.00 1.01
## [505,] 1.00 1.01
## [506,] 1.00 1.02
## [507,] 1.00 1.02
## [508,] 1.00 1.02
## [509,] 1.00 1.02
## [510,] 1.00 1.02
## [511,] 1.00 1.02
## [512,] 1.00 1.02
## [513,] 1.00 1.02
## [514,] 1.00 1.02
## [515,] 1.00 1.02
## [516,] 1.00 1.02
## [517,] 1.00 1.02
## [518,] 1.00 1.02
## [519,] 1.00 1.02
## [520,] 1.00 1.02
## [521,] 1.00 1.02
## [522,] 1.00 1.02
## [523,] 1.00 1.01
## [524,] 1.00 1.01
## [525,] 1.00 1.01
## [526,] 1.00 1.01
## [527,] 1.00 1.01
## [528,] 1.00 1.01
## [529,] 1.00 1.01
## [530,] 1.00 1.01
## [531,] 1.00 1.01
## [532,] 1.00 1.01
## [533,] 1.00 1.01
## [534,] 1.00 1.01
## [535,] 1.00 1.01
## [536,] 1.00 1.01
## [537,] 1.00 1.02
## [538,] 1.00 1.02
## [539,] 1.00 1.02
## [540,] 1.00 1.01
## [541,] 1.00 1.01
## [542,] 1.00 1.01
## [543,] 1.00 1.01
## [544,] 1.00 1.01
## [545,] 1.00 1.01
## [546,] 1.00 1.01
## [547,] 1.00 1.01
## [548,] 1.00 1.01
## [549,] 1.00 1.01
## [550,] 1.00 1.01
## [551,] 1.00 1.01
## [552,] 1.00 1.01
## [553,] 1.00 1.01
## [554,] 1.00 1.01
## [555,] 1.00 1.02
## [556,] 1.01 1.03
## [557,] 1.01 1.03
## [558,] 1.01 1.03
## [559,] 1.01 1.02
## [560,] 1.00 1.02
## [561,] 1.00 1.01
## [562,] 1.00 1.01
## [563,] 1.00 1.01
## [564,] 1.00 1.00
## [565,] 1.00 1.00
## [566,] 1.00 1.00
## [567,] 1.00 1.00
## [568,] 1.00 1.00
## [569,] 1.00 1.00
## [570,] 1.00 1.00
## [571,] 1.00 1.00
## [572,] 1.00 1.00
## [573,] 1.00 1.00
## [574,] 1.00 1.00
## [575,] 1.00 1.00
## [576,] 1.00 1.00
## [577,] 1.00 1.00
## [578,] 1.00 1.00
## [579,] 1.00 1.00
## [580,] 1.00 1.01
## [581,] 1.00 1.01
## [582,] 1.00 1.01
## [583,] 1.00 1.01
## [584,] 1.01 1.03
## [585,] 1.00 1.02
## [586,] 1.00 1.01
## [587,] 1.00 1.00
## [588,] 1.00 1.00
## [589,] 1.00 1.00
## [590,] 1.00 1.00
## [591,] 1.00 1.00
## [592,] 1.00 1.00
## [593,] 1.00 1.00
## [594,] 1.00 1.00
## [595,] 1.00 1.00
## [596,] 1.00 1.01
## [597,] 1.01 1.02
## [598,] 1.01 1.02
## [599,] 1.01 1.02
## [600,] 1.01 1.02
## [601,] 1.01 1.02
## [602,] 1.01 1.02
## [603,] 1.00 1.02
## [604,] 1.00 1.01
## [605,] 1.00 1.01
## [606,] 1.00 1.01
## [607,] 1.00 1.01
## [608,] 1.00 1.01
## [609,] 1.00 1.00
## [610,] 1.00 1.00
## [611,] 1.00 1.00
## [612,] 1.00 1.01
## [613,] 1.00 1.01
## [614,] 1.00 1.00
## [615,] 1.00 1.00
## [616,] 1.00 1.00
## [617,] 1.00 1.00
## [618,] 1.00 1.00
## [619,] 1.00 1.00
## [620,] 1.00 1.01
## [621,] 1.00 1.01
## [622,] 1.00 1.00
## [623,] 1.00 1.01
## [624,] 1.00 1.01
## [625,] 1.00 1.01
## [626,] 1.00 1.01
## [627,] 1.00 1.01
## [628,] 1.00 1.00
## [629,] 1.00 1.00
## [630,] 1.00 1.01
## [631,] 1.00 1.00
## [632,] 1.00 1.01
## [633,] 1.00 1.01
## [634,] 1.00 1.01
## [635,] 1.00 1.01
## [636,] 1.01 1.01
## [637,] 1.00 1.01
## [638,] 1.00 1.01
## [639,] 1.00 1.01
## [640,] 1.00 1.01
## [641,] 1.00 1.01
## [642,] 1.00 1.00
## [643,] 1.01 1.01
## [644,] 1.01 1.01
## [645,] 1.01 1.01
## [646,] 1.01 1.01
## [647,] 1.01 1.02
## [648,] 1.01 1.02
## [649,] 1.01 1.02
## [650,] 1.01 1.01
## [651,] 1.01 1.01
## [652,] 1.01 1.01
## [653,] 1.01 1.01
## [654,] 1.00 1.00
## [655,] 1.00 1.00
## [656,] 1.00 1.00
## [657,] 1.00 1.01
## [658,] 1.00 1.01
## [659,] 1.00 1.01
## [660,] 1.01 1.02
## [661,] 1.01 1.02
## [662,] 1.00 1.01
## [663,] 1.00 1.01
## [664,] 1.00 1.01
## [665,] 1.00 1.01
## [666,] 1.00 1.00
## [667,] 1.00 1.00
## [668,] 1.00 1.00
## [669,] 1.00 1.00
## [670,] 1.00 1.00
## [671,] 1.00 1.00
## [672,] 1.00 1.00
## [673,] 1.00 1.00
## [674,] 1.00 1.00
## [675,] 1.00 1.01
## [676,] 1.01 1.01
## [677,] 1.01 1.02
## [678,] 1.01 1.02
## [679,] 1.01 1.02
## [680,] 1.01 1.02
## [681,] 1.01 1.02
## [682,] 1.01 1.02
## [683,] 1.01 1.02
## [684,] 1.01 1.02
## [685,] 1.01 1.03
## [686,] 1.01 1.02
## [687,] 1.01 1.01
## [688,] 1.00 1.00
## [689,] 1.00 1.01
## [690,] 1.00 1.01
## [691,] 1.00 1.01
## [692,] 1.00 1.01
## [693,] 1.00 1.01
## [694,] 1.00 1.01
## [695,] 1.00 1.00
## [696,] 1.00 1.00
## [697,] 1.00 1.00
## [698,] 1.00 1.00
## [699,] 1.00 1.00
## [700,] 1.00 1.01
## [701,] 1.00 1.01
## [702,] 1.01 1.01
## [703,] 1.01 1.01
## [704,] 1.01 1.01
## [705,] 1.01 1.01
## [706,] 1.00 1.00
## [707,] 1.00 1.01
## [708,] 1.00 1.01
## [709,] 1.01 1.02
## [710,] 1.00 1.02
## [711,] 1.00 1.02
## [712,] 1.00 1.02
## [713,] 1.00 1.01
## [714,] 1.00 1.01
## [715,] 1.00 1.01
## [716,] 1.00 1.00
## [717,] 1.00 1.00
## [718,] 1.01 1.01
## [719,] 1.00 1.00
## [720,] 1.00 1.00
## [721,] 1.00 1.01
## [722,] 1.01 1.01
## [723,] 1.01 1.01
## [724,] 1.01 1.01
## [725,] 1.01 1.01
## [726,] 1.01 1.02
## [727,] 1.02 1.05
## [728,] 1.03 1.06
## [729,] 1.03 1.07
## [730,] 1.02 1.06
## [731,] 1.03 1.06
## [732,] 1.03 1.08
## [733,] 1.03 1.10
## [734,] 1.03 1.09
## [735,] 1.03 1.09
## [736,] 1.04 1.11
## [737,] 1.04 1.10
## [738,] 1.03 1.06
## [739,] 1.02 1.04
## [740,] 1.01 1.02
## [741,] 1.01 1.01
## [742,] 1.01 1.03
## [743,] 1.01 1.02
## [744,] 1.01 1.01
## [745,] 1.01 1.01
## [746,] 1.00 1.01
## [747,] 1.00 1.01
## [748,] 1.01 1.03
## [749,] 1.01 1.05
## [750,] 1.01 1.05
## [751,] 1.01 1.04
## [752,] 1.01 1.04
## [753,] 1.01 1.03
## [754,] 1.01 1.03
## [755,] 1.01 1.03
## [756,] 1.01 1.03
## [757,] 1.00 1.02
## [758,] 1.00 1.02
## [759,] 1.00 1.01
## [760,] 1.00 1.01
##
## Multivariate psrf
##
## 1.37
gelman.plot(lista)
# Tentar descobrir como colocar label no gelman plot